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We present a formalism for investigating the interaction between p-mode oscillations and 
convection by analyzing realistic, three-dimensional simulations of the near-surface layers of the 
solar convection zone. By choosing suitable definitions for fluctuations and averages, we obtain 
a separation that retains exact equations. The equations for the horizontal averages contain 
one part that corresponds directly to the wave equations for a 1-D medium, plus additional 
terms that arise from the averaging and correspond to the turbulent pressure gradient in the 
momentum equation and the divergence of the convective and kinetic energy fluxes in the internal 
energy equation. These terms cannot be evaluated in closed form, but they may be measured in 
numerical simulations. The additional terms may cause the mode frequencies to shift, relative to 
what would be obtained if only the terms corresponding to a 1-D medium were retained — most 
straightforwardly by changing the mean stratification, and more subtly by changing the effective 
compressibility of the medium. In the presence of time dependent convection, the additional 
terms also have a stochastic time dependence, that acts as a source of random excitation of the 
coherent modes. In the present paper, we derive an expression for the excitation power and test 
it by applying it to a numerical experiment of sufficient duration for the excited modes to be 
spectrally resolved. 
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ABSTRACT 



The near-surface layers of the Sun are of crucial 
importance for the properties of the solar p-mode 
oscillations (see the recent conference proceedings; 
Brown 1993; Hoeksema et al. 1995; Ulrich et al. 
1995; Antia & Chitre 1996; Pijpers et al. 1997). 
The upper turning points of the p-modes are lo- 
cated in these layers and this is where the modes 
are excited and damped. This is also where the so- 
lar convection zone gives way to the visible solar 



The thin superadiabatic layer at the top of the 
solar convection zone is characterized by large fluc- 
tuations in the thermodynamic variables. As il- 
lustrated by detailed numerical simulations of the 
solar surface layers (Nordlund 1982, 1985; Steffen 
et al. 1989; Nordlund & Dravins 1990; Nordlund 
& Stein 1991a; Steffen & Freytag 1991; Stein & 
Nordlund 1989, 1991, 1994, 1998; Atroshchenko & 
Gadun 1994; Solanki et al. 1996), the fluctuation 
amplitudes peak just below the visible surface. 



where the temperature ranges from 5000 to 10,500 
K, and the logarithmic fluctuations of the density 
and pressure, A In p and A In P, are of the order of 
unity. The velocity amplitudes are large through- 
out the photosphere, with rms Mach numbers of 
the order of 0.3 and peak Mach numbers exceeding 
unity in a small fraction of the volume (Nordlund 
& Stein 1991b; Stein & Nordlund 1998). 

The layers with large amplitude fluctuations (at 
recent meetings referred to as the 'muck region') 
may be expected to influence the solar p-mode os- 
cillations in several ways. First, these are the lay- 
ers wlicirc! most of the random excitation of modes 
is expected to occur (Stein 1967, 1968; Goldreich 
& Keeley 1977; Goldreich & Kumar 1988, 1990; 
Stein & Nordhmd 1991; Bogdan et al. 1993; Gol- 
dreich et al. 1994; Musielak et al. 1994). Second, 
the average vertical stratiflcation of this region 
cannot be assumed to be in hydrostatic equilib- 
rium; the motions will contribute an additional 
'turbulent pressure', that adds to the normal gas 
pressure and hence tends to elevate the surface 
layers. The large amplitude, non-linear fluctua- 
tions makes even the deflnition of appropriate av- 
erage values a non-trivial exercise. Further, due to 
the extreme temperature sensitivity of the opacity 
in the surface layers the emergent solar luminos- 
ity is produced by an average state that diflfers 
noticeably from that of a corresponding one di- 
mensional model. Third, because of the presence 
of the fluctuations, the wave propagation prop- 
erties of the medium will in general be different 
than for a homogeneous medium. We follow Balm- 
forth (1992a) and refer to the effects due to mean 
structure changes as 'extrinsic' (or 'model') effects, 
and those that are caused by changes in the wave 
propagation properties of the medium as 'intrinsic' 
('modal' or 'mode physics') effects. 

A problem with reproducing the frequencies of 
the modes that have upper turning points in these 
layers has indeed been known since the early days 
of helioseismology (Christensen-Dalsgaard et al. 
1996b; Christensen-Dalsgaard 1988; Christensen- 
Dalsgaard et al. 1988). Improvements in the calcu- 
lation of the equation of state (Mihalas et al. 1990; 
Rogers et al. 1996) did not improve the situation, 
but instead rather sharpened the significance of 
the discrepancy between the observed and calcu- 
lated oscillation frequencies. 

The discrepancy between the observed and the- 



oretical mode frequencies is primarily a function 
of frequency and is nearly independent of degree, 
£, for small £. It is small for the lowest frequen- 
cies, and grows to significant values for frequen- 
cies approaching the cut-off frequency of the solar 
photosphere (approximately 5 mHz). This shows 
that the cause of the discrepancy resides in layers 
to which the low frequency modes hardly pene- 
trate, but where the high frequency modes have 
a significant amplitude. Thus, the source of the 
discrepancy must lie near the solar surface, in the 
outer layers of the solar cavity. 

Stocrhastic excitation of p-modes has been 
demonstrated in numerical simulations of convec- 
tion in the solar surface layers (Stein et al. 1988; 
Stcffcn 1988; Stein et al. 1989; Stein & Nordlund 
1991; Bogdan et al. 1993). In principle, the vari- 
ous contributions to stochastic excitation, damp- 
ing and frequency shift may be directly measured 
in such numerical simulations. Alternatively, one 
may instead extract information about the model 
structure and mode propagation properties from 
the numerical simulations, and carry that infor- 
mation over to standard envelope and mode calcu- 
lation procedures. The advantage with the latter 
method is that one is not limited to the studying 
the sparse spectrum of modes that are excited in 
a small box. The main purpose of the present 
paper is to present a formalism for interpreting 
quantities used in such 1-D calculations as suit- 
able averages of quantities that may be measured 
in 3-D numerical simulations. 

Stein & Nordlund (1991) used a simplified ver- 
sion of the formalism presented here in an ini- 
tial study of mode excitation and Rosenthal ct al. 
(1999); Rosenthal et al. (1998) also used a simpli- 
fied version to analyze the effect of 1-D/3-D model 
differences on the frequencies of radial modes, 
showing that on the one hand the model differ- 
ences may account for the majority of the fre- 
quency discrepancy, but that on the other hand 
the modal effects also are significant. 

In the present paper we develop a formalism 
for analyzing the interaction of convection with 
purely radial oscillations by choosing an exact de- 
composition into horizontal averages and fluctua- 
tions. Since most p- modes are nearly radial near 
the solar surface, the analysis covers the lower or- 
der behavior of non-radial modes as well. 

In Section 2 we present the separation of vari- 
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ables that wc have chosen to work with. In Sec- 
tion 3 we use this formahsm to analyze how the 
interaction of convection with the oscillations can 
cause mode excitation, mode damping and fre- 
quency shifts, and in Section 4 we show by ex- 
plicit application to a numerical experiment that 
the expression for the mode excitation produces 
estimates of mode power that are consistent with 
what is actually observed in the experiment. 

In order to verify the formula for the mode ex- 
citation it is necessary to use a numerical exper- 
iment of sufficiently long duration for the excited 
modes to be spectrally resolved, because; this al- 
lows the mode excitation power to be estimated di- 
rectly from a measurement of the mode power and 
the mode lifetime (obtainable from the fuU-width- 
at-half-maximum of the mode energy). This ne- 
cessitates a relatively low spatial resolution, which 
precludes a direct comparison with the solar mode 
excitation power (currently well established, see 
for example Roca Cortes et al. (1999)). 

In a subsequent paper (Paper II) we make use 
of simulations with higher spatial resolution (and 
correspondingly shorter duration), that allow us to 
make detailed comparisons with helioseismic data, 
and to explore details of the processes that domi- 
nate the stochastic mode excitation. 

2. Formalism 

Ideally, wc would like to split the fluid equations 
into one set describing the oscillations and one set 
describing convection. By the very nature of the 
problem, such a separation cannot be complete; 
if the convection is to excite the oscillations, it 
must give rise to a source term in what would oth- 
erwise be equations describing ideal, radial wave 
motions. And if convection is to have an effect 
on the (complex) frequencies of the modes, there 
must be a possibility for convection to affect the 
compressibility of the gas; in other words to have a 
coherent response with a phase lag that in general 
may be expected to vary with height. In addition, 
the convection is able to affect the frequencies in a 
more trivial manner, namely by changing the over- 
all stratification relative to whatever reference 1-D 
model one happens to compare with. 

We are thus content with, and indeed looking 
for, a separation consisting of wave-like equations 
with additional terms that are related to the pres- 



ence of convection. Such a separation is possibly 
not unique, but below we present one possibility. 

One factor influencing our choice of separation 
was that we are here concerned with a case that 
is, in a sense, opposite to the more common case 
of small amplitude 3-D fluctuations on top of large 
means. We have large amplitude 3-D fluctuations 
due to convection on top of a mean with small 
amplitude coherent p-mode fluctuations. Thus we 
prefer to work with the actual equations, rather 
than some truncated expansion. We retain explic- 
itly those terms that cannot be worked out ana- 
lytically, and subsequently measure them in the 
numerical simulations. In doing so, we wish to 
make the split such that these terms have a well 
defined physical meaning, and are numerically well 
conditioned. 

2.1. Notation and definitions 

To achieve the goals set out above, it is cru- 
cial to choose an appropriate set of definitions for 
fluctuations and averages. We use the convention 
that per-unit- volume quantities arc written in up- 
per case and per-unit-mass quantities are written 
in lower case. 

For the present case of purely radial motions, 
we define straight horizontal averages F and corre- 
sponding residuals F for per-unit-volume variables 
F, 

FM = {F)^y (1) 
F{x,y,z,t) = F{x,y,z,t)-F{z,t), (2) 

and density weighted horizontal averages / with 
residuals / for per-unit-mass variables /, 

= {Pf)j{p).y (3) 

f{x,y,z,t) = fix,y,z,t)^ f{z,t). (4) 

Thus, horizontal averages vanish for residuals of 
per-unit-volume variables, 

{F)^y = 0, (5) 

while for the residuals of per-unit-mass variables 
it is the mass weighted horizontal averages 

{Pf),y = (6) 

that vanish. 
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Because of the large number of occurrences of 
horizontal averages, we henceforth drop their ex- 
plicit xy subscript (O^.^^ ()) and retain sub- 
scripts only on time averages and ensemble 
averages (O^^J. 

As a consequence of Eq. 3, products of two and 
three per-unit-mass variables obey 

ipfg) = {pm + PfS (7) 

{pfgh) = {pfgh) + pfgh + 

{pgh)f+{pfh)g+{pfg)h (8) 

2.2. 1-D averaged equations 

To derive a set of hydrodynamic equations for 
the horizontal averages, we apply the formalism of 
the previous subsection to the three-dimensional 
equations of mass, momentum and energy conser- 
vation, 



The terms on the RHS may be interpreted as the 
kinetic energy flux of the radial motions, the ad- 
vection of convective kinetic energy density by the 
radial motions, the Pu flux associated with the 
turbulent pressure (cf. below), and the radial ki- 
netic energy flux associated with the horizontal 
fluctuations. 

An equation for conservation of total internal 
plus kinetic energy, e = Ci + Ck, may be obtained 
by adding the time derivative of the kinetic energy 



density, ek 

ei, 

d 



1„,2 



to that for the internal energy 



(pe) = -V • (peu + Pu + Frad + Fvisc) + pu • g, 

(15) 



dt 

where F^isc is the viscous flux 



j 



(16) 



-V-(pu), (9) 
-V ■ (puu -a) -VP- pg, (10) 
-V- (pciu) -P(V-u) 

-V- Frad + Qdiss , (H) 



where ei is the internal energy, P = P{p, Ci) is the 
gas pressure, and F^ad is the radiative energy flux. 
Qdiss is the viscous dissipation 



Horizontal averaging of the above equations 
yields equations that depend on height and time, 
but not on the horizontal coordinates: 



d_ 
dt 



g-^ipU.) = 



iP) 



d_ 



{pUz), 



(17) 



d ,-2 



•Pg - cTzz) + gp 



(18) 



Qdiss — ^ ^ ^ij ^ij ) 



(12) 



where is the symmetric part of the strain tensor 
dui/dxj, and aij is the viscous stress tensor, cry = 

pVSij . 

By Eqs. 7-8, the total kinetic energy {^pu^) 
splits into the kinetic energy of the radial motion 
and of the horizontal fluctuations, 



{^pu ) = ^pu^ + {-pu ) 



(13) 



and the radial kinetic energy flux {^pv? Uz) splits 
into 



{\pu'^ Uz) 



1 — 2- , / -i "2\ - 
= 2^zUz + {^PU )Uz 

+ {piiDuz + {\pU^Uz). (14) 



dt 



(pe) 



[ (-^conv -^rad) 

az 

+ Odiss + (ii-VP)], (19) 



d ___ 
- ^ {pe-Uz + Pu^ ) g/ow^ 
az 

(-^conv ~t~ -^rad 



dz 

-^visc)] 



(20) 



conv are defined 



where^p = (p), pu^_ = {pu^), Pg 
Pg + Pt, and pe = (pe) . Pt and 
below (Eqs. 21-22). 

Equations 17-20 correspond closely to the hy- 
drodynamic equations for a stratified, homoge- 
neous 1-D medium. The terms that are not brack- 
eted are exactly what appear in the 1-D equations. 
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whereas the extra, bracketed terms are due to the 
convective motions. These are the gradient of the 
turbulent pressure 



the divergence of the convective flux 



((pei + Pg)u,), 



(21) 



(22) 



and the divergence of the kinetic energy flux asso- 
ciated with convection, 



(23) 



Here Uz is the vertical velocity relative to a frame 
of reference moving with velocity Uz = {uzp)/{p); 
we refer to this frame of reference as 'pseudo- 
lagrangian'. Note that, in general, {uz) ^ 0. 

The 

-Prad term should, in principle, not be 
bracketed, since it also may appear in a 1-D 
medium. It is, however, so intimately related 
to the other flux terms that we prefer to place 
it inside the brackets. Furthermore, because of 
the strong non-linearities involved, the value of 
Frad niay differ substantially between 1-D and 3-D 
models with the same average structure. 
A stationary state that has Uz =Q obeys 



d 

-{P^ + P,-azz) 



9P 



(24) 



dz 



{Fconv + Jkin + -FVisc + i^rad) = (25) 



Note that the horizontally averaged pressure is 
in general not the same as the pressure for the av- 
erage density and energy of which it is a function. 

Because of the correlation of the fluctuations and 
the contribution from higher order terms, 



g — Pg{p, e) -|- ^ {pe) + higher order terms 
" ' (26) 



dpde 

is in general not equal to Pg{p, e). 



2.3. Pseudo-Lagrangian 1-D Equations 

For the purpose of studying adiabatic (or near 
adiabatic) wave mode fluctuations, it is more rele- 
vant to write the horizontally averaged equations 
in the pseudo-lagrangian reference frame (moving 
with velocity Uz)- The operator 



D d d 



(27) 



picks up the time variation in that frame of refer- 
ence. For a per-unit-mass variable / that satisfies 



|(p/) = -V.F, 



(28) 



the change to a pseudo-lagrangian frame results in 



^(p/) = -W-F^Uz^ipf) 



d 



= _V.(F-p/u)-/)/— (n,)(29) 

i.e., the subtraction of the average flux of pf from 

within the divergence operator, and the addition 
of a term of the form —pf-^{uz), that accounts 
for the dilation or concentration of the (per unit 
volume) quantity p f , corresponding to stretching 
or compression of the coordinate system. 

In terms of the -^() operator the averaged 1-D 
equations become 



D 

m 



(P) 



D 
Dt 



{pUz) 



D 
Dt 



(pe) 



-p§;i^^) (30) 

-^(^g - ^zz) + 9p-pUz-^{Uz) 

4(A)], (31) 

d , __ __ 9 ,_ , 
(Puz ) + gpuz - pe—(Uzj 



dz 
d 



dz 



[ q^{Fqot[iy ~f~ -^kin ~t~ -^visc ~(~ -^rad)]- 

(32) 

The transformation between per- unit- volume and 
per-unit-mass pseudo-lagrangian time derivatives 
is 



^(P7)-/^(P-) 

The dilation term that is present in the per-unit- 
volume formulation thus drops out in the per-unit- 
mass formulation, because a fixed amount of mass 
is under consideration — the mass element pdz is 
indeed suitable for depth integration of pseudo- 
lagrangian per-unit-mass quantities. 

In per-unit-mass variables, the equation of mo- 
tion and the energy equations for horizontal aver- 
ages thus become 

P;^(^.) = -^(^g + A-^z.) + 5P,(34) 
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d 



[ o (-^conv ^ -^rad) 

+ Qdiss + (u-VP}], (35) 
d -_ 

d - 

[ — -Q^ (-Fconv + -Frad 
+ fkin + ^;isc)]. (36) 



In what follows, we use the traditional notation 
5Pg, Sp, etc., to distinguish pseudo-lagrangian per- 
turbations from Eulerian ones. 

3. Interactions between oscillations and 
convection 

We now use this formaUsm to calculate the work 
done on radial oscillation modes by convection. 
We then test the resulting expression on the ver- 
tical resonant modes excited in a numerical simu- 
lation of solar convection. 

In the presence of a coherent mode, the time 
variation of the additional, convcctivc terms in 
Eqs. 17-20 is partly in unison with the mode, 
reflecting the coherent response of convection to 
the presence of the mode. The coherent response 
may again be divided into one part that is in 
phase with the mode (appearing as the real part 
of a Fourier transform), and one part that is in 
quadrature with the mode (the imaginary part of 
a Fourier transform). The imaginary part of the 
coherent component causes an exponential damp- 
ing (or growth) of a trapped mode, and the real 
part causes a frequency shift. In analogy with 
simpler situations we also refer to these parts as 
"adiabatic" and "non-adiabatic" . There is also an 
incoherent contribution, corresponding to the ran- 
dom variation of the convection, that contributes 
to these averages, and produces stochastic mode 
excitation and damping. 

The equation describing the time evolution of 
the mode kinetic energy may be obtained from Eq. 
34: 



d - - 



+ (Pg + Pt - a,^) 



dz 



Integrating over time and depth, we obtain 
J dM^ul = Jdtjdz ((5Pg + SPt 



-5a, 



. du. 



dz 

+ j dt j dz puzg 

~t- [ ••• ] boundaries* (3^) 

where, [ ... ] boundaries dcuotcs the boundary contri- 
butions that result from integrating the divergence 
form of the equation. For suitably chosen bound- 
ary conditions, these contributions vanish. Specif- 
ically, the displacement may be chosen to vanish 
at the lower boundary, and the pressure may be 
chosen sufficiently small at the \ipper boundary. 
The work done by gravity vanishes if there are no 
net mass displacements in the model. 

The Jdtjdz {SPg + 5P^ - S(f,,_)du,/dz part of 
the work integral represents the PdV work done 
by the gas pressure, the turbulent pressure, and 
the mean viscous stress. (The mean viscous stress 
may be expected to be negligibly small in stellar 
atmospheres and envelopes, and should be small 
also in numerical simulations.) The duz/dz fac- 
tor is equal to —Dlnp/Dt = D\nV/Dt, where 
V is the speciflc volume. Pressure perturbations 
must be out of phase with density perturbations 
in order to contribute to the work integral. In a 
diagram showing SP against 5p this corresponds 
to open curves, with a counter clockwise sense of 
orientation. 

The signs are as to be expected from basic phys- 
ical principles; the mode kinetic energy increases if 
the pressure is larger during expansion than dur- 
ing compression. 

3.1. Coherent and incoherent fluctuations 

The PdV work integral, Eq. 38, is 

W = Jdtjdz SP-^. (39) 

Here 6P is the pseudo-lagrangian total pressure 
fluctuation (neglecting the small viscous stress 
contribution) 



SP = SP^ + 5Pt. 



(40) 



(37) 
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and ^ is the displacement, which is related to the 
velocity by 



negligible. The (a) x (1) contribution, for example. 



^ = J dt'u,, 
and to the density variations by 
Dlnp di 



Dt 



dz' 



(41) 



(42) 



The displacement and the pressure fluctuations 
can be split into two parts: one the coherent, sinu- 
soidal, modal part and the other the incoherent, 
random part: 



5P = SP^ + 5Pr. 



(43) 
(44) 



The pressure fluctuations can be further split into 
an adiabatic part, proportional to the density fluc- 
tuations, 

6lnP'^=Ti{z)5lnp = -T,{z)^, (45) 
and a non-adiabatic part, the remainder: 

(5A; = 



dz 



Pr^iz)^+SPj\ (46) 



SPr = -Pri(2)^+^P/, (47) 

so that, 

5P = JP„" + SP^" + SPr" + SP". (48) 

Thus, the work integral can be expanded into the 
product of 4 pressure fluctuation terms multiplied 
by 2 displacement terms, 

W = j j dtdz 

[SP^" + dPj- + SPr" + (5P/] 



dz 



dz 



dtdz 



[(a) + (6) + (c) + (d)][(l) + (2)]. (49) 

Consider the contribution of each possible pair to 
the work. The integral as a whole may be expected 
to display a "random walk" behavior; i.e., its ex- 
pectation value grows with total time of integra- 
tion. Contributions that are bounded are therefore 



(a) X (1) cx 



djoj dju 

dz dz 



D 

m 
1 

2 



di^ 
dz 



di^ 
dz 



t2 



tl 



is a total integral, so it is negligible. 

(c) X (2) ( 



dir 

dz 



d£,r dir 


dz 


dz 


D 






Dt 


[K 


1 




^d^ 


2 




^ dz 



(50) 
(51) 

(52) 

(53) 
(54) 

(55) 



tl 



so its contribution is also negligible, 
(a) X (2) + (c) X (1) oc 







+ pri 


oz 


dz 






^rr# 

Dt 


( di^ 


di^ 


) 


V dz 


dz 




^PFi 


' di^ 


dir' 


t2 


dz 


dz 


tl 



djr dj^ 
dz dz 



(56) 



and also gives a negligible contribution. Note that 
this step depends on using the same Fi in the def- 
initions of the adiabatic modal and random parts 
of the pressure fluctuation (Eqs. 46 and 47). 



(6) X (1) oc 5P^ 



d^ 

dz 



(57) 



represents the linear driving/ damping of the 
mode. It is the balance of this term with the 
stochastic driving that determines the amplitude 
of the mode. 



(6)x(2) oc 5P."^ 
oz 



(58) 



oc P(51nP>r^i<51nP/, (59) 

is a stochastic driving term. However, it is small 
in comparison to the next term, because it is pro- 
portional to the non-adiabatic coherent pressure 
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fluctuation which is small since the mode ampli- 
tude growth time is many periods. 



(c) X (1) cx SPr 



dz 



(60) 



oc PSlnPj'ujY^^SlnPl, (61) 

is the dominant stochastic driving term. It is pro- 
portional to the non-adiabatic random pressure 
fluctuations which arc large and the adiabatic co- 
herent pressure fluctuations which are also large, 
since they provide the mode restoring force. 



(d) X (2) oc 5Pr 



dz' 



(62) 



does not represent any work on the mode (it does 
not contain any factor representing the mode). 

3.2. Stochastic excitation 

Using the results from the previous subsec- 
tion, wc may derive expressions that measure 
the stochastic excitation and the linear damping 
in numerical simulations. After proper scaling, 
such measurements yield estimates of global ex- 
citation power and damping that may be com- 
pared directly with estimates from observations 
Roca Cortes et al. (1999), and with estimates from 
analytical theories (e.g., Balmforth 1992b; Goldre- 
ich et al. 1994). 

The mode energy per unit surface area (at r = 
R) is 



(63) 



where M^^ is (by definition) the mode mass, and 
Voj = ^ui{R) is the mode velocity amplitude at the 
reference radius R. The change in a mode's kinetic 
energy over a time interval At is 

Ai;^ = f dt [ dr 6P ^ 

J At Jr or 



where 



= [ dtES W^{t) 

J At 

W^{t) = J^dr dP 



(64) 
(65) 



is the instantaneous work integral for the given 
mode's displacement, normalized with the square 
root of the mode's energy, and 

C{iv,At) = [ dt W^{t) (66) 

J At 

is the work integrated over the time interval At. 

For small changes of amplitude \{j Vui+AVui. 
Equation 63 with the definitions in Eq. 64 give, 

AV^ = ^M~^C{uj,At). (67) 

For a particular realization of the stochastic driv- 
ing, there is a complex phase factor e*'^ between 
A^ and V^^. The ensemble average of the mode 
energy E^^ + AE^^ at t + At is proportional to the 
ensemble average of -|- A"(4,|^ over all phases. 
In this averaging the linear term (K<j^K;)ens 
ishes, and one obtains from the quadratic term 



A(|K.^|)ens = (|AK,|\ 



(68) 



\-M^-^C{oj,At)f) ,(69) 

^ ens 



SO that 



At 



^u.A(|yJ|)_ 



At 



1 



^^^JC(a;,At) 



(70) 
• (71) 



E£ C(w, At) = At Re[ / dr 6P*iuj- 

J r 

= uAt Im[ / dr 5P* 



In terms of Fourier transforms over the time inter- 
val At, the integrated work 

^1 
dr ^ 

3* d^ui 1 

dr^ 

= ujAtlm[W^] (72) 

is proportional to the Fourier amplitude of the pro- 
jection of the stochastic pressure fluctuations onto 

The expectation value of the square of the imag- 
inary part is half of the power of this random func- 
tion in the frequency interval Av = 1/At, since 

(Mens = (I MW] |2)ens + (I MW] l'),^,, and 

the ensemble averages for the real and imaginary 
parts are equal for a function with randomly dis- 
tributed phases, so 



(73) 
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Prom Eqs. 70, 72 and 73, we obtain 



4. Examples 



e 

|W^j')e„s 



8E^ ' 



fdrSP:^\^) (74) 

•/r ens 



and with the definition Eq. 63 of the mode energy, 
the expectation value of the stochastic excitation 
power is 



At 4Aiyldr\U^p{^)^- 
From Eq. 57, the damping rate is 

1 dE^ _ IdrSP^ 



at ^^2j^ar\u'p (i) 



2 ' 



or, in terms of Fourier transforms, 

1 dE^ _2ldr lm[5P* ^] 



E^ dt 



(75) 



(76) 



(77) 



Equation 75 expresses the excitation power per 
unit surface area in terms of the power density of 
the stochastic pressure fluctuations, SP^;. If these 
are measured in a numerical simulation covering 
a small {L x L) patch of the solar surface, the re- 
sult must be properly scaled in order to estimate 
global excitation power. Taking the results at face 
value would correspond to assuming a periodic, co- 
herent repetition of the pressure fluctuations over 
the entire solar surface area AwR"^. If the actual 
pressure fluctuations are instead uncorrelated over 
scales larger than the scale of the numerical model, 
the global spectral power density is smaller by a 
factor of Li^ / {AttB?). The integrated excitation 
power is then obtained by multiplying Eq. 75 with 
the horizontal surface area of the numerical model, 
i^, rather than by the total surface area, A-kB?. 
(Note that this does not imply that the scale L 
enters into the final result, since the expectation 
value of the fluctuations of the mean pressure 5Pu, 
drops correspondingly, if the size of the numeri- 
cal model is increased beyond the size over which 
pressure fluctuations are correlated.) 



To demonstrate that the proposed formalism 
is useful, both from a theoretical and a practical 
point of view, we provide two example applica- 
tions; a theoretical discussion of deviations from 
hydrostatic equilibrium, and a practical applica- 
tion to a numerical data set. 

4.1. Deviations from hydrostatic equilib- 
rium 

The equations for the horizontal averages (ei- 
ther Eqs. 17-20 or the equivalent Lagrangian equa- 
tions 30-32) do not assume anything about the 
horizontally averaged (barred) quantities, other 
than that they are instantaneous, horizontal av- 
erages. 

Assuming now that a stationary reference state 
exists, with (Wz)^ = 0, one may ask what happens 
if initially the barred variables deviate from the 
equilibrium state. The initial state may then be 
taken as an initial condition for a time integration 
of Eqs. 30-32. 

Ignoring in a first analysis the lagrangian time 
variation of the terms that represent the convec- 
tivc flux, the radiative flux, and the turbulent 
pressure, one has a set of equations that are di- 
rectly analogous to the ordinary, one dimensional 
hydrodynamic equations, and admits exactly the 
same type of solution; the initial non-equilibrium 
state slumps towards the equilibrium state, over- 
shoots, and a scries of oscillations ensues. If the 
boundaries are closed, and there are no dissipative 
terms, the oscillations will continue undamped, 
and consist of a superposition of the eigenmodes 
of the system. 

If the boundaries allow wave energy to escape, 
and / or there arc other dissipative effects, the 
eigenmodes will be damped, and the solution will 
eventually settle to the equilibrium solution. 

Likewise, if one considers the actual, time- 
dependent lagrangian response of the convcctive 
and radiative fluxes, and of the turbulent pressure, 
these will in general also contribute to the damp- 
ing (or possibly self-excitation) of the eigenmodes 
of the system. 

In conclusion, an initial state that is not in 

hydrostatic equilibrium may be thought of as a 
superposition of (in general damped) eigenmodes. 
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and such a system will typically perform damped 
oscillations around a stationary reference state. 

All of this might be analyzed by making use 
of the formalism put forward here, at any chosen 
level of realism; one may replace the bracketed 
terms that are not available in closed form with 
theoretical estimates, or one may retain them as 
they are, and evaluate them from numerical sim- 
ulations. In fact, we have performed numerical 
experiments of just this type by subjecting snap 
shots from 3-D simulations to large amplitude per- 
turbations of the vertical equilibrium, in order to 
study in particular the damping properties of the 
resulting set of oscillations. Results of the analysis 
will be published in a forthcoming paper. 

4.2. Excitation rates and energy losses in 
a numerical experiment 

As a direct test of the expression for the excita- 
tion rate, Eq. 74, we have evaluated the expression 
numerically, using data from a numerical experi- 
ment covering 43 hours of solar time at a resolution 
of 63 X 63 X 63 additional details about this ex- 
periment are given elsewhere. The length of this 
run is sufficient for a direct measurement of mode 
line widths, and hence for a direct determination 
of the right hand side of the expression 



mode 



ode 



(78) 



where the mode energy life time Tmode is related 
to the full width at half maximum of the mode 
energy spectrum Az/pwHM by 



27rrmodeA!/FWHM = 1- 



(79) 



In Fig. 1 we compare estimates of the excitation 
power based on Eq. 75 with the actual energy loss 
for the three radial modes that are observed in the 
experiment. The energy loss of each mode was 
computed from Eq. 78 by fitting Lorentzians to 
the excess mode power above the background noise 
(cf. Paper II, Fig. 1), using Smodc and Tmodc as pa- 
rameters in the fit. The excitation power was eval- 
uated from Eq. 75, with mode displacement fac- 
tors and ^tj, obtained by using Fourier trans- 
forms to project out the three modes (cf. Fig. 2). 
Note that the actual amplitudes of these modes do 
not enter into the estimate of the excitation power, 
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Fig. 1. — The excitation power (small pluses) in 
the neighborhood of the frequencies of the three 
radial modes that are excited to measurable ampli- 
tudes in the box. The average excitation is shown 
as diamonds, and the average energy loss per unit 
time is shown as squares. The results are from a 
numerical experiment with realistic opacities and 
equation of state, covering 43 hours of solar time 
at a resolution of 63 x 63 x 63 (Gcorgobiani et al. 
2000). Note that, because of the relatively low 
resolution required for such a long run to be af- 
fordable, the results are not directly comparable 
to solar values. The mode mass of the box also 
differs significantly from the solar mode mass at 
the same frequency. 
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Fig. 2. — The velocity mode shapes for the three 
radial oscillations, at 2.5, 3.9, and 5.6 mHz, that 
are observed in the 43 hour simulation. For com- 
parison, the corresponding modes from Model S of 
Christensen-Dalsgaard et al. (1996a), normalized 
to the same mode mass within the box, are also 
shown (plus symbols) for the two modes whose fre- 
quencies are below the acoustic cut-off frequency. 




Fig. 3. — The power-spectrum of the non- 
adiabatic pressure fluctuations, Eq. 47, in a layer 
close to the surface of the 43 hour simulation. 



because of the normalization with the mode ki- 
netic energy in Eq. 75. The non-adiabatic pres- 
sure fluctuations at the mode frequencies, 5P*, 
were obtained from the Fourier transform of the 
non-adiabatic pressure fluctuations (Fig. 3), aver- 
aging over a neighborhood of the mode frequen- 
cies to obtain the values indicated with diamonds 
in Fig. 1. For further details on how the formulae 
are evaluated see the appendix of Paper II. 

Note that the spectrum of non-adiabatic pres- 
sure follows a power law and shows no particular 
features at the resonant mode frequencies. This 
illustrates that the non-adiabatic fluctuations are 
mostly incoherent. Conversely, since the modes 
are only weakly damped, the coherent fluctuations 
at the mode frequencies are mostly adiabatic. At 
non-mode frequencies, the fluctuations are a mix 
of adiabatic and non-adiabatic, incoherent fluctu- 
ations (of which only the non-adiabatic ones con- 
tribute stochastic work). 

5. Concluding remarks 

We have shown how to decompose the fluid 

variables into horizontal averages and fluctuations, 
and how to use this decomposition to separate the 
equations for the radial p-modes (Eqs. 30-32) from 
the full equations (Eqs. 9-11). The equations for 
the radial modes are similar to those for a 1-D 
stratifled medium with the addition of two extra 
terms ~ the gradient of the turbulent pressure and 
the gradient of the convective plus kinetic energy 
fluxes. There are two additional subtle differences 
from the 1-D equations - the radiative flux may 
differ significantly from that of a 1-D model with 
the same mean structure and the gas pressure is 
not the same as the pressure for the average den- 
sity and internal energy. 

We have used the decomposed fluid equations 
to derive an expression for the stochastic excita- 
tion rate of the radial p-modes (Eq. 75) in terms 
of the PdV work of the non-adiabatic, incoherent, 
random, convcctively produced pressure fluctua- 
tions and the mode compression, and tested the 
expression by applying it to a numerical experi- 
ment of sufficient length to measure the mode en- 
ergy loss directly from the observed mode power 
and spectral line width. 

The price that had to be paid for obtaining 
a sufficiently long duration experiment (43 solar 
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hours in this case) was that the niiiiK^rical res- 
olution could not be very large. From previous 
convergence investigations we know that various 
aspects of the numerical models depend to a quite 
varying degree on the numerical resolution (Stein 
& Nordlund 1998). The thermal mean structure, 
for example, is very robust, while the peak in rel- 
ative turbulent pressure near the surface depends 
more sensitively on numerical resolution. 

In the subsequent paper (Stein & Nordlund 
2000, Paper II), we use higher resolution models 
to determine the mode excitation power more ac- 
curately, and compare it directly with helioseismic 
results. The duration of these experiments are not 
sufficient to allow the excited modes to be resolved 
in frequency, but this is not necessary to determine 
the mode excitation power from Eq. 75. In addi- 
tion, the high resolution experiments are used to 
reveal the spatial locations that contribute most of 
the excitation power, and to investigate the nature 
of the mechanism responsible for the excitation. 

Using the type of formalism put forward in 

the present paper it is also possible to analyze 
the mode physics (intrinsic) contributions to mode 
damping and frequency shifts, but this is work for 
the future. 
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Figure Captions 



Fig. 1. The excitation power (small pluses) in 
the neighborhood of the frequencies of the three 
radial modes that are excited to measurable ampli- 
tudes in the box. The average excitation is shown 
as diamonds, and the average energy loss per unit 
time is shown as squares. The results arc from a 
numerical experiment with realistic opacities and 
equation of state, covering 43 hours of solar time 
at a resolution of 63 x 63 x 63 (Georgobiani et al. 
2000). Note that, because of the relatively low 
resolution required for such a long run to be af- 
fordable, the results are not directly comparable 
to solar values. The mode mass of the box also 
differs significantly from the solar mode mass at 
the same frequency. 

Fig. 2.- The velocity mode shapes for the three 
radial oscillations, at 2.5, 3.9, and 5.6 mHz, that 
are observed in the 43 hour simulation. For com- 
parison, the corresponding modes from Model S of 
Christensen-Dalsgaard et al. (1996a), normalized 
to the same mode mass within the box, are also 
shown (plus symbols) for the two modes whose fre- 
quencies are below the acoustic cut-off frequency. 

Fig. 3. — The power-spectrum of the non- 
adiabatic pressure fluctuations, Eq. 47, in a layer 
close to the surface of the 43 hour simulation. 
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